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We study theoretically the erosion threshold of a granular bed forced by a viscous fluid. We first 
introduce a novel model of interacting particles driven on a rough substrate. It predicts a continuous 
transition at some threshold forcing 0 C , beyond which the particle current grows linearly J ~ 8 — 0 C , 
in agreement with experiments. The stationary state is reached after a transient time t conv which 
diverges near the transition as t conv ~ \0 — 0 C \~ z with 2 « 2.5. The model also makes quantitative 
testable predictions for the drainage pattern: the distribution P(a) of local current is found to be 
extremely broad with P(a) ~ J/a, spatial correlations for the current are negligible in the direction 
transverse to forcing, but long-range parallel to it. We explain some of these features using a scaling 
argument and a mean-field approximation that builds an analogy with g-models. We discuss the 
relationship between our erosion model and models for the depinning transition of vortex lattices in 
dirty superconductors, where our results may also apply. 


Erosion shapes Earth’s landscape, and occurs when 
a fluid exerts a sufficient shear stress on a sedimented 
layer. It is controlled by the dimensionless Shields num¬ 
ber 0 = E/(p p — p)gd , where d and p p are the particle 
diameter and density, and p and E are the fluid density 
and the shear stress. Sustained sediment transport can 
take place above some critical value 9 C [1-3], in the vicin¬ 
ity of which motion is localized on a thin layer of order 
of the particle size, while deeper particles are static or 
very slowly creeping [4-6]. This situation is relevant in 
gravel rivers, where erosion occurs until the fluid stress 
approaches threshold [7]. In that case, predicting the flux 
J of particles as a function of 0 is difficult, both for tur¬ 
bulent and laminar flows [4, 8]. We focus on the latter, 
where experiments show that: (i) in a stationary state, 
J oc (6 — 6 C )P with /3 « 1 [4, 6, 9, 10], although other ex¬ 
ponents are sometimes reported [3], (ii) transient effects 
occur on a time scale that appears to diverge as 9 —> 9 C 
[4, 6] and (iii) as 9 —> 9 C the number of moving particles 
vanishes, but not their characteristic speed [4, 10]. 

Although a continuous description of erosion appears 
successful for 9 9 C [5, 9, 11], it should not apply for 9 —► 
9 C . In the latter regime, an erosion/deposition model was 
proposed in [4], where one assumes that a 0-dependent 
fraction of initially mobile particles evolve over a frozen 
static background, which contain holes. In this view, 9 C 
occurs when the number of holes matches the number of 
initially moving particles. This phenomenological model, 
which assumes no interactions between mobile particles, 
captures (i,ii,iii) qualitatively well. This success is sur¬ 
prising: due to the frozen background, one expects mo¬ 
bile particles to take favored paths and to eventually 
clump together into ’’rivers”, thus avoiding most of the 
holes. Models including this effect as well as particle in¬ 
teractions [12, 13] have been introduced in the context 
of the depinning transition of vortex lattice in dirty su¬ 
perconductors. They lead to a sharp transition for the 


flux at some finite forcing, but /3 ss 1.5. Moreover, there 
are currently no predictions for the spatial organization 
of the flux near threshold, although this property is in¬ 
dicative of the underlying physics, and could be accessed 
experimentally. 

In this letter we introduce a model of interacting parti¬ 
cles forced along one direction on a disordered substrate. 
Particle interactions based on mechanical considerations 
are incorporated. Such model recovers (i,ii,iii) with /3 = 1 
and an equilibration time t conv ~ \9—9 c \~ 2 5 . In addition, 
we find that (a) the spatial distribution of local flux a is 
extremely broad, and follows P(a) ~ 1/a and (b) spatial 
correlations of flux are short-range and very small in the 
lateral direction, but are power-law in the mean flow di¬ 
rection. We derive f} = 1 and explain why P(a) is broad 
using a mean-field description of our model, leading to 
an analogy with q-models [14, 15] used to study force 
propagation in granular packings. 

Model: we consider a density n of particles on a frozen 
background, n should be chosen to be of order one, but 
its exact value does not affect our conclusions. The back¬ 
ground is modeled via a square lattice, whose diagonal in¬ 
dicates the direction of forcing, referred to as “downhill”. 
The lattice is bi-periodic, of dimension L x W, where L 
is the length in downhill direction and W the transverse 
width. Each node i of the lattice is ascribed a height 
hi £ [0,1], chosen randomly with a uniform distribution. 
Lattice bonds i —> j are directed in the downhill direc¬ 
tion, and characterized by an inclination = hi — hj. 
We denote by 9 the amplitude of the forcing. For an 
isolated particle on site i, motion will occur along the 
steepest of the two outlets (downhill bonds) [16], if it 
satisfies 9 + 9,^j > 0. Otherwise, the particle is trapped. 

However, if particles are adjacent, interaction takes 
place. First, particles cannot overlap, so they will only 
move toward unoccupied sites. Moreover, particles can 
push particles below them, potentially un-trapping these 
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or affecting their direction of motion. To model these 
effects, we introduce scalar forces fi->j on each outlet of 
occupied sites, which satisfy: 

fi—tj — max(fj/- ¥ i T @i—*j + 9, 0) (1) 

where f is the force on the input bond j' —> i along 
the same direction as i —> j , as depicted in Fig. 1. Eq.(l) 
captures that forces are positive for repulsive particles, 
and that particle i exerts a larger force on toward site j 
if the bond inclination is large, or if other particles 
above i are pushing it in that direction. From our analysis 
below, we expect that the details of the interactions are 
not relevant, as long as the direction of motion of one 
particle can depend on the presence of particles above it- 
an ingredient not present in [12, 13]. 



FIG. 1. Illustration of the model. Small circles indicate lattice 
sites, particles are represented by discs in yellow, or green if 
motion occurred between t (left) and t + 1 (right). The black 
arrow is in the downhill direction. Solid lines indicate outlet 
with positive forces. If a particle has two outlets with positive 
forces, the larger (smaller) one is colored in red (blue). 

We update the position of the particles as follows, see 
Fig. 1 for illustration. We first compute all the forces in 
the system. Next we consider one row of W sites, and 
consider the motion of its particles. Priority is set by con¬ 
sidering first outlets with the largest fi^j and unoccu¬ 
pied downhill site j. Once all possible moves ( > 0, 

j empty) have been made, forces are computed again in 
the system, and the next uphill row of particles is up¬ 
dated. When the L rows forming the periodic system 
have all been updated, time increases by one. 

For given parameters 9 , n we prepare the system via 
two protocols. In the “quenched” protocol, one consid¬ 
ers a given frozen background, and launch the numerics 
with a large 0 and randomly placed particles - parame¬ 
ters are such that the system is well within the flowing 
phase. Next, 9 is lowered slowly so that stationarity is al¬ 


ways achieved. We also consider the “Equilibrated” pro¬ 
tocol: for any 9 , particles initial positions are random. 
Dynamical properties are measured after the memory of 
the random initial condition is lost. We find that us¬ 
ing different protocols does not change critical exponents, 
but that the quenched protocol appears to converge more 
slowly with system size. Below we present most of our 
results obtained from the “equilibrated” protocol with 
W = AVL [13], and n = 0.25 unless specified. 



FIG. 2. Average current J versus 9 — 9 c in log-log scale for the 
(a) “equilibrated” and (b) “quenched” protocols, for which 
9 C = 0.164 ± 0.002 and 9 C = 0.172 ± 0.002 respectively- a 
difference plausibly due to finite size effects. The black solid 
lines with slope one indicate the linear relation J oc 9 — 8 C . 
(c) Density of conducting sites p s versus 9 — 9 C for the “equi¬ 
librated” protocol, (d) p a curves collapsed by rescaling 9 — 9 C 
with L 1 /", where v = 3.0 ± 0.2. 


Results: Once the steady state is reached, we measure 
the average current of particles J and the number density 
of sites carrying a finite current p s . Measurements of 
both quantities indicate a sharp dynamical transition at 
some 9 C below which J = 0 and p s = 0 as L —> oo, see 
Fig. 1. 9 C can be accurately extracted by considering the 
crossing point of the curves p s {9) as L is varied, yielding 
9 C = 0.164 ± 0.002 for the equilibrated protocol. In the 
limit L —> oo our data extrapolates to: 

J{9) ~9-9 c for 9 > 9 C (2) 

P,(0)=e(p-9 C ), (3) 

where O is the Heaviside function. Eq.(2) corresponds to 
/3 = 1, whereas Eq.(3) indicates that all sites are visited 
by particles in the flowing phase. Introducing the expo¬ 
nent p s {9) ~ {9 — 9c) 1 , this corresponds to 7 = 0. The 
collapse of Fig. 2(d) shows how convergence to Eq.(3) 
takes place as L —> 00 , from which a finite size scaling 
length £ ~ (9 — 9 C )~ 1J with v ss 3 can be extracted. 
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Criticality is also observed in the transient time t conv 
needed for the current to reach its stationary value. Fig. 3 
reports that t conv ~ \9 — 0 C \~ z with 2 « 2.5 on both 
sides of the transition. A similar exponent was observed 
numerically in [17]. 



FIG. 3. Transient time tconv v.s. 9. For a given realization, 
tconv is defined as the smallest time for which J(f) — J < 
y/Var{J) where Var(J) = liinT->oo ^ — J) 2 . The 

gray dashed lines correspond to t conv ~ 1 9 — 9 C \~ 2 ' 5 . 

The spatial organization of the current in steady state 
can be studied by considering the time-averaged local 
current cr, on site i. or the time-averaged outlet current 
The spatial average of each quantity is J. Fig. 4 
shows an example of drainage pattern, i.e. one realization 
of the map of the ■ 


FIG. 4. Examples of drainage pattern just below 9 C (Left) 
and above (Right). The black arrow shows the downhill direc¬ 
tion. The thickness of the lines represents in logarithmic 
scale. A few examples showing splitting events are magnified 
on the left. Here W = 45 and L = 128, and J > 0 even below 
9 C due to finite size effects. 

To quantify such patterns, we compute in Fig. 5(a) the 
distribution P(cr) of the local current cr; for various mean 
current J. We observed that: 

( 4 ) 


where r w 1 and / is a cut-off function, expected since 
in our model cq < 1. Eq.(4) indicates that P(cr) is re¬ 
markably broad. In fact, the divergence at small <7 is so 
pronounced that a cut-off <r m i n must be present in Eq.(4) 
to guarantee a proper normalization of the distribution 
P(cr), although we cannot detect it numerically. 



FIG. 5. Distribution of the site current P(ff) in steady state 
for given average currents J of (a) the erosion model (L = 256, 
W = 64) and (b) our mean-field model (IF = 1600). 

Next, we compute the spatial correlation of the local 
current in the transverse direction C't(x), defined as: 

C T (x) = ((c TiOi+x) - J 2 )/(Wr) ~ J 2 ) (5) 

where the site i and i+x are on the same row, but at a dis¬ 
tance x of each other. Here the brackets denote the spa¬ 
tial average, whereas the overline indicates averaging over 
the quenched randomness (the hi s). Fig. 6(a) shows that 
no transverse correlations exist for distances larger that 
one site. However, long-range, power-law correlations are 
observed in the longitudinal direction, as can be seen by 
defining a longitudinal correlation function Cl(v), where 
y is the vertical distance between two sites belonging to 
the same column. We find that Cl(v) ~ 1/v^ #c, but 

that Cl(v) decays somewhat faster deeper in the flowing 
phase, as shown in Fig. 6(b). 



FIG. 6. (a) Transverse current correlations Ct at 9 C and (b) 
longitudinal current correlation Cl at 9 C and at 9 — 9 C = 0.25 
for L = 256 (dashed line). 

A scaling relation: we now derive a relationship be¬ 
tween the exponents /3 characterizing J and 7 charac- 




P(a) = Jcr~ T f(a) 
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terizing p s . It holds true for both protocols, but is pre¬ 
sented here in the “quenched” case. Near threshold, at 
any instant of time the density of moving particles is 
J <C n < 1, thus most of the particles are trapped and 
will move only when a mobile particle passes by. As 0 
is decreased by some SO , a finite density of new traps 
8m ~ SO is created. If these traps appear on the region 
of size p s where mobile particles flow, they will reduce 
the fraction of mobile particle by SJ = p s Sm ~ PsSO, 
which implies: 


= 7+1 


( 6 ) 


Eq. (6) shows that the result 0 = 1 is a direct consequence 
of the fact that in our model, all sites are explored by mo¬ 
bile particles for 0 > 0 C , a result which is not obvious. 
In the dirty superconductor models of [13, 18], this is 
not the case and for the “equilibrated” protocol 0 > 1 
was found. We argue that this difference comes from the 
dynamical rules chosen in [13, 18], according to which 
“rivers” forming the drainage pattern never split: their 
current grows in amplitude in the downhill direction, un¬ 
til it reaches unity. In these models the drainage pattern 
thus consists of rivers of unit current, avoiding each other, 
and separated by a typical distance of order 1 / J. Our 
model behaves completely differently because rivers can 
split, as emphasized in Fig. 4. This comes about because 
the direction taken by a particle can depend on the pres¬ 
ence of a particle right above it, as illustrated in case A 
of Fig. 1. This effect is expected to occur in the ero¬ 
sion problem due to hydrodynamic interactions or direct 
contact between particles, and may also be relevant for 
superconductors. 

Mean-field model: we now seek to quantify the effect of 
splitting. Its relevance is not obvious a priori, as splitting 
stems from particle interactions, and may thus become 
less important as the fraction of moving particles van¬ 
ishes as J —> 0. To model this effect we consider that 
the current cq on a site i is decomposed in its two out¬ 
lets as (Ji = qcri + (1 — q)(ii, where q is a random vari¬ 
able of distribution r)(q). If there were no splitting then 
r](q) = hS(q) + 1 — q). Here instead, we assume that 

77 (g) = \5{q — J) + ^<5(1 — J — q). This choice captures 
that the probability of splitting is increased if more mov¬ 
ing particles are present, and can occur for example if 
two particles flow behind each other, as exemplified in 
case A of Fig. 1. Next, we make the mean field assump¬ 
tion that two adjacent sites i and j on the same row are 
uncorrelated, P(ai,aj) = P{ai)P{aj). We then obtain 
the self-consistent equation that P(cr) must be equal to: 


J dq 1 dq 2 dcr 1 d(72r)(qi)r](q2)P(<Ji)P(cr2)5(q 1 <T 1 + q 2 a 2 - q) 

(7) 

This mean-field model belongs to the class of g-models in¬ 
troduced to study force propagation [14, 15]. It is easy to 
simulate, and some aspects of the solution can be com¬ 
puted. Numerical results are shown in Fig. 5(b). The 


result obtained for P{a) is very similar to Eq.(4) that de¬ 
scribes our erosion model: P(cr) is found to be power-law 
distributed (although r = 3/2 instead of r = 1) where 
with an upper cutoff at <r max ~ 1, and P(a) oc J. 

These results are of interest, as they explain why P(a) 
is very broad, and is not dominated by sites displaying no 
current at all (which would correspond to a delta function 
at zero) even as J —> 0, thus confirming that 7 = 0. 
They can be explained by taking the Laplace transform 
P of Eq.(7). One then obtains a non-linear differential 
equation for P, from which it can be argued generically 
that r = 3/2 [15]. We have performed a Taylor expansion 
of P around zero, which leads to relationship between the 
different moments of the distribution P(ct). From it, we 
can show that P(cr) oc J and cr ma x ~ 1. We also find that 
the cut-off of the divergence of P(cr) at small argument 
follows (7 m in ~ J 1 /( T_1 ). 

Conclusion: we have introduced a novel model for over¬ 
damped interacting particles driven on a disordered sub¬ 
strate. It predicts a dynamical phase transition at some 
threshold forcing 0 C , and makes quantitative predictions 
for various quantities including the particle current and 
the drainage pattern. The latter could be tested exper¬ 
imentally in erosion experiments [3, 4, 6, 9] by tracking 
particles on the surface [4] to reconstruct the spatial or¬ 
ganization of current. Another interesting set-up are col¬ 
loids at an interface, pinned by a random environment 
generated by a rough charged surface [19]. Numerics 
support the existence of a dynamical transition in this 
system where flow localizes on channels [20], which may 
fall in the universality class of our model. 

Note that our model assumes that particles are over¬ 
damped, and that their inertia is negligible. We expect 
inertia to lead to hysteresis and make the transition first 
order, as observed on inertial granular flows down an 
inclined plane [21], although this effect may be small in 
practice, as supported by experiments [22]. We did not 
consider non-laminar flows, nor temperature (that can 
be relevant for colloids). Both effects should smooth the 
transition, and lead to creep even below 0 C . 

Finally, it has been proposed that the erosion thresh¬ 
old is a dynamical transition very similar to the jamming 
transition that occurs when a bulk amorphous material 
is sheared [6]. If our model holds, this is not the case: 
due to the presence of the free interface, long-range elas¬ 
tic interactions between mobile particles are absent. In 
recent theoretical descriptions of the jamming transition 
such interactions are central both for soft [23] and hard 
particles [24, 25]. 
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